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We address the folding induced by differential growth in soft layered solids via an elementary model 
that consists of a soft growing neo-Hookean elastic layer adhered to a deep elastic substrate. As 
the layer/substrate modulus ratio is varied from above unity towards zero we find a first transition 
from supercritical smooth folding followed by cusping of the valleys to direct subcritical cusped 
folding, then another to supercritical cusped folding. Beyond threshold the high amplitude fold 
spacing converges to about four layer thicknesses for many modulus ratios. In three dimensions the 
instability gives rise to a wide variety of morphologies, including almost degenerate zigzag and triple¬ 
junction patterns that can coexist when the layer and substrate are of comparable softness. Our 
study unifies these results providing understanding for the complex and diverse fold morphologies 
found in biology, including the zigzag precursors to intestinal villi, and disordered zigzags and triple¬ 
junctions in mammalian cortex. 


I. INTRODUCTION 

Some biological tissues have evolved into folded mor¬ 
phologies to accommodate large surfaces in small vol¬ 
umes; folding of the mammalian brain is associated with 
expansion of its outer cerebral cortex, and the folded ep¬ 
ithelium in the gut provides large surface area for ab¬ 
sorption of nutrients. Pattern formation in morphogene¬ 
sis is typically attributed to chemical instabilities within 
Turing’s react ion-diffusion paradigm [T], but increasing 
experimental evidence shows that these folded morpholo¬ 
gies actually emerge from smooth precursors via mechan¬ 
ical instabilities driven by differential growth of layered 
tissues [HE]. Recent studies have considered mechanical 
descriptions of brain folding ms], epithelial folding in 
the gut 019 ] and carcinomas m, fingerprint formation 
m, and circumferential buckling in tubular organs and 
tumors [HHH]. 

We consider the most elementary model for soft tissue 
folding: a growing elastic layer with shear modulus /ii 
and relaxed thickness h adhered to an infinite substrate 
with shear modulus /i2 = Ui/O- The layer’s growth in¬ 
duces compression, ultimately causing it to buckle Hl- 
[29] . If 77 ^ 1 the layer buckles supercritically into smooth 
wrinkled patterns p!5ll2Q] with wavelength I ^ ^ h 

explicable via small strains and linear stability analysis 
[21]. If 77 ^ 1, the substrate does not deform, but af¬ 
ter sufficient growth the layer’s free surface invaginates 
into cusped furrows called sulci or creases via an exotic 
supercritical non-linear instability that eludes analytic 
treatment |3Qll36] . However, most cases of biological in¬ 
terest involve a growing soft tissue on a comparably soft 
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substrate - a mechanically intriguing regime [22H29] that 
must span a transition from smooth to cusped folding. 
Previous theoretical studies have focused on two dimen¬ 
sions and have considered post-bifurcation behavior when 
approaching this transition from the stiff-layer smooth 
wrinkling side |26ll28] , and the linear stability of circular 
layers either side of the transition [29]. Here we use lin¬ 
ear analysis and numerics to explore folding around the 
transition itself, including beyond threshold behavior in 
two and three dimensions. We discover the transition 
occurs via an unexpected region of subcritical cusped 
folding and is accompanied by an unusually rich set of 
three-dimensional pattern morphologies, which we argue 
underpin the mechanics and diversity of invaginations in 
developmental biology. 


II. RESULTS 


We examine the stability of a layer occupying 0 < 2: < 
a atop a substrate occupying z < 0. The layer has un¬ 
dergone growth, swelling, or pre-strain such that, with¬ 
out adhesion, it would undergo a relaxing deformation 
G = didig{gx^ gy, gz)^ becoming a flat slab with thick¬ 
ness h = a/gz. We model both tissues as incompress¬ 
ible neo-Hookean elastic solids so that, if this reference 
state is subject to a displacement u and corresponding 
deformation gradient F = / + Vu, the elastic energy 
per unit volume is |/i2Tr (F • F^) in the substrate and 
|/iiTr(F.G-i-G-^-F^) in the layer, while incom¬ 
pressibility requires Det(F) = Det(G) = 1. 
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A. Linear instability 


zero wavelength instability happens continuously at 


We first consider uniaxial growth, G = diag(^, 1, 1 /^), 
and consider the layer’s linear stability to infinites¬ 
imal isochoric perturbations, u = f{z) sm{kx)z + 
{f'{z)/k)cos{kx)ii. The full calculation is in Appendix 
[B] here we restrict ourselves to an overview. The lin¬ 
earized Euler-Lagrange equations for the solids require 
^ 4 j ( 4 ) _ ^ 1 )/^^) + k'^f = 0 (with ^ = 1 in the 

substrate), solved by f{z) = a4e^^ -h + a2e^^/^ + 

in the layer and f{z) = e^^{bi + zkb2) in the 
substrate. Imposing material continuity at the interface 
and natural boundary conditions at the free surface and 
interface we conclude the layer becomes linearly unstable 
when g satisfies 


/ be fees" \ 
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where b = 1 g^. Fig. plots this threshold as a func¬ 

tion of wavelength for different modulus ratios. Minimiz¬ 
ing this threshold over wavelengths gives the wavelength 
and threshold of the first unstable mode as a function 
of 77, plotted in Fig. [^-c. As expected the instabil¬ 
ity reproduces the long-wavelength wrinkling limit for 
77 1 and the zero wavelength Biot instability [22] at 

9 = 9 b = 1 . 839 ... for 77 ^ 1 . The transition from finite to 



FIG. 1. Elastic linear instability of a growing slab on a 
substrate, showing a transition from finite to zero wavelength 
instability at modulus ratio gt ~ 0.525. (a) Growth g required 
for instability at wavelength I relative to a deformed thickness 
a for a range of moduli ratios, (b-c) Wavelength and threshold 
growth of the first unstable mode as a function of 77 . Dashed 
lines indicate the asymptotic results for stiff layers. 
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However, the zero wavelength Biot instability is never ob¬ 
served: it is always preceded by a point of nonlinear insta¬ 
bility to infinitesimal cusped sulci [ 32 ] . In incompressible 
neo-Hookean materials this occurs at a uniaxial stretch 
As ~ 0.647 [ 32 ], corresponding to ^s = 1 /As ~ 1 . 55 . 
Thus we expect finite wavelength wrinkling as predicted 
by linear stability for 77 > 77s ~ 1 . 14 , where the linear 
threshold intersects gg (see Fig. [^), and direct sulcifica- 
tion at gs for 77 < 77s. 


B. Folding in two dimensions 

We verify these predictions and explore folding be¬ 
yond threshold with numerical simulations (Appendix 
0- starting with uniaxial growth in two dimensions. 
When g > gs the surface indeed buckles smoothly at 
the linear instability threshold. This buckling amplifies 
surface compression in the valleys, producing the sulcifi- 
cation strain there shortly thereafter. The valleys then 
sulcify in a secondary bifurcation which occurs before gs 
up to 77 « 10 (Fig. Cusps form first on the surface 
and then on the interface between the layer and substrate 
at high g (Fig. [^). However, the interface in real sys¬ 
tems may not be strictly sharp, suppressing interfacial 
cusping. 

When 77 < 77s sulcification sets in at gs^ before the lin¬ 
ear instability. For 77 near gs the energy of the system at 
g = gs'i^ minimized by a cusped fold that has a finite am¬ 
plitude (see the layer with 77 = 1 and ^ = 1.55 ~ ^s in Fig. 
HI), implying that sulcification is a nonlinear subcritical 
instability when 77 is near 775. The instability is known 
to be supercritical in both the stiff layer wrinkling limit 
77 ^ 1 and in the 77 —> 0 pure sulcifications limit [ 33 l | 36 ] , 
so this region of sub-criticality between the two limits is a 
surprise, stemming from hybridization of the two modes. 
In Fig. [^we show hysteresis loops for the instability at 
three values of 77 characterizing the three regimes: linear 
and super-critical for 77 ^ 1, non-linear and subcritical 
for g ^ 1 and non-linear and supercritical for g < 0 . 5 . 
Analogous subcritical instabilities were recently identi¬ 
fied in two-dimensional elastic systems with jointly com¬ 
pressed layer and substrate nil EH]. These systems differ 
from ours particularly strongly in the region of subcriti- 
cality: indeed at 77 = 1 the distinction between layer and 
substrate vanishes in the jointly compressed system, so 
the wavelength scales with the total depth of the system, 
whereas in our system the layer is distinguished by its 
growth even at 77 = 1 and the wavelength is always finite 
despite the substrate being infinite. 

We optimized the wavelengths of the folds for mini¬ 
mum energy (Fig.|^) finding that in the smooth-folding 
regime, 77 > 77^, the wavelength I decreases past threshold 
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FIG. 2. Folding of a uniaxially growing layer in two dimensions, (a) Morphologies as a function of modulus ratio ij and growth 
ratio g = Qx. The diagram includes the linear instability curve (blue line, dotted for 77 < r]s ~ 1-14 where the sulcification 
threshold is met before the linear instability) and numerically determined thresholds for cusp formation (black line). (b) Energy 
minimizing wavelengths as a function of g for a range of moduli ratios 77. Hollow and solid circles correspond to smooth and 
cusped folds, respectively. The blue curve is the wavelength of linear instability (the dotted end is above the sulcification 
threshold gs ~ 1.55). (c) Simulated minimum energy states near threshold (left column) and far from threshold (middle and 
right column) for 77 = 0.5, 1 and 4. Color indicates isotropic stress: dark red < — 5/i, dark blue > 5/i (/x = gi in the layer 
and /X = in the substrate), (d) Cortical folds in a ferret brain (braincatalogue.org). (e) Cross-sections of previllous folds in 
intestinal epithelium of a chick embryo at days 13 and 14 [40] (credits: A. E. Shyer). 


(in Eulerian coordinates), which is similar to folding of 
stiff layers m- In contrast, when f] < rjs the optimal 
wavelength increases past threshold. Consequently the 
optimal wavelength settles to about four times the un¬ 
deformed layer thickness over a wide range of 77 when g 
becomes large. This agrees with the folding patterns in 
the brain and the gut I3IH1, and is a hallmark of mechan¬ 
ical compressive folding. 

The low 77 folds are notably asymmetric, with a deeply 
cusped upper surface and a wavy bottom surface, in con¬ 
trast to nearly top-bottom symmetric folds at 77 > 77^ 
that take a sawtooth-like form at high amplitudes (Fig. 
[^). This agrees again with the observations in folded 
tissues; the brain has a low modulus ratio between the 
cortex and sublayers (77 ^ 1 [2]) leading to thick and 
deeply cusped asymmetric folds (Fig. [^) [ 39 ], whereas 
narrower folds with a sawtooth profile ^ig. |^) are seen 
during gut development, owing to a relatively stiff ep¬ 


ithelium that folds on a softer mesenchyme [3 So]. We 
note, however, that in a developing chick gut the epithe¬ 
lium and mesenchyme are both confined by outer layers 
of smooth muscle, which contributes to folding at early 
stages [ 3 ]- 


C. Folding in three dimensions 

If G = (lmg{gx^ {gxdy)) is not uniaxial, but we 
nevertheless consider x — z plane-strain deformations 
then, as seen in eqn. |B 3 [ the elastic energy depends 
on two dimensionless parameters, rj/g‘^ and rjg^gy. Our 
plane-strain calculations thus hold with the identifica¬ 
tions 77 ^ gyf] and g gx-s/^- In particular, for trans¬ 
versely isotropic growth, g = gx = gyi Ihe sulcification 
threshold is at gx^/^ — g^^^ ~ 1 . 55 , i.e, gs ~ 1 . 34 , 
and the linear threshold intersects gs at rjs ~ 0 . 86 . This 
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FIG. 3. Hysteresis loops for three different modulus ratios, demonstrating the three different folding regimes, (a) When 
?7 = 0.5 sulcification sets in supercritically at p ~ 1.55, producing a non-linear supercritical instability, reminiscent of pure 
sulcihcation. (b) When 77 = 1 sulcification occurs subcritically ai g = gs, immediately producing a deep cusped fold, while in 
unfolding the fold also disappears sub-critically at a lower value of g. Folding is non-linear and subcritical. (c) When 77 = 4 
folding is linear and supercritical, reminiscent of pure wrinkling. 


f] = 4 





plane-strain calculation captures the full linear stability 
of the system since the Fourier modes that would make 
up a truly three-dimensional pattern cannot interact at 
the level of linear stability. 

Before examining folded states in three dimensions, we 
recall previous results in the stiff and soft layer limits: In 
stiff layers (77 ^ r]s) a hexagonally ordered pattern of 
dents forms near threshold [191 120 ] (although analyses 
based on the plate theory predict a checkerboard pat¬ 
tern in the stiff film limit [T8l[T9]) and a zigzag (herring¬ 
bone) pattern is the preferred far-from-threshold pattern 
naEo]. In contrast, in clamped soft layers (77 = 0 ) we 
previously identified a square pattern of line-like sulci 
with alternating orientations as the optimal near thresh¬ 
old pattern and a hexagonal array of triple-junction sulci 
as the optimal far-from-threshold pattern [ 36 ]. There 
should thus be transitions between the two sets of pat¬ 
terns as the layer stiffness is reduced from large toward 
zero. 

Whereas folding of relatively stiff layers can be stud¬ 
ied even at large overstresses using experimental setups 
based on swelling elastomers mun], folding of soft lay¬ 
ers is harder to realize in experimental model systems be¬ 
cause large swelling ratios are needed to reach far-from- 
threshold states. Patterns of line-like sulci and triple¬ 
junction sulci have been produced experimentally in soft 
layers adhered to rigid substrates pusniisi], but very 
few experiments have explored three dimensional folding 
patterns in the regime of modest stiffness ratios 0 ( 24 ]. 

We simulate three dimensional folding patterns in¬ 
duced by transversely isotropic growth of the layer on 
periodic domains that contain one unit cell of the pat¬ 
tern with energetically optimized wavelengths (Appendix 
A). We use external guiding forces to inititate the pat¬ 
terns that, once formed, are stable in the neighborhood 
of 77,5f-space where they are the energy minimizing pat¬ 
terns. Optimal near-threshold patterns, including a rect¬ 
angular array of line-like sulci at 77 ^ 77^ and smooth 


hexagonal dents at 77 > 77^, are shown on the left of 
Fig. [^. These patterns minimize the energy only at 
e/ec close to unity, where we define the relative growth 
strain e/ec = conveniently quantify distance from 

the threshold (e = g — ^ and €c = gc — ^ where gc is 
the growth ratio at the threshold). As in two dimen¬ 
sions, folding is subcritical when 77 is near 77^, evidenced 
by a finite amplitude pattern at the threshold g = gs 
when f] = rjs (Fig. &)• For far-from-threshold states 
(e/ec ^ 2) we construct zigzag patterns as well as pat¬ 
terns of hexagonal triple-junctions, discovering that these 
two patterns are almost equally good energy minimizers 
(Fig. [^). However, the triple-junctions are slightly bet¬ 
ter at low 77 and/or high e/ec, while zigzags become rela¬ 
tively better toward large 77, as expected. Fig. shows 
a phase diagram summarizing these findings. 

The optimal wavelengths of the zigzags and triple¬ 
junctions are shown in Fig.[^ and[^, respectively. As in 
the uniaxial case, the wavelengths decrease with g when 
T] > rjs 1 but remain nearly steady when g ^ rjs 1 result¬ 
ing in an approximate collapse of the wavelengths toward 
large g such that the width of the folds is again roughly 
four times the undeformed layer thickness. This means 
that further expansion of the layer is accomodated by 
increasing amplitude. The optimal x, ^-aspect ratio of 
zigzags decreases with decreasing 77 such that the zigzags 
in soft layers are relatively shorter than their stiff-layer 
counterparts m- 

Given the nearly equal energies of the zigzag and triple¬ 
junction folds, it is not surprising that systems with mod¬ 
est modulus ratios, such as the folded cerebral cortex 
(Fig. 1^) or its gel layer mimics 0 , show disordered mix¬ 
tures of zigzags and triple junctions (in contrast, epithe¬ 
lium of the gut shows ordered zigzags, see Fig. [J, owing 
to a higher modulus ratio and anisotropy 0 ). By simu¬ 
lating unaided folding of the biaxially growing layer in a 
large simulation domain whose width is multiple times 
the wavelength of the patterns (Appendix A), we see 
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FIG. 4. Three dimensional folding of an equibiaxially growing layer, (a) Energy minimizing periodic patterns near threshold 
(left) and far from the threshold (right, zigzag and triple-junction patterns are shown for each 77). (b) Diagram of energy 
minimizing ordered patterns as a function of 77 and g. The symbols (Z’s for zigzags, Y’s for triple-junctions, o’s for smooth 
dents, etc.) indicate the minimum energy pattern at each simulated point. The blue curve is the threshold for linear instability 
or sulcification. The dashed lines are contours for growth strains e/cc. (c) Deformation energies of zigzags and triple-junctions 
relative to the energies of flat reference states. (d,e) Energetically optimal half-wavelengths [lx, ly, and Ih are indicated in 
(a)] of zigzags and triple-junctions, respectively. Thick black curves indicate the wavelength of linear-instability, (f) Zigzags 
in epithelium of a developing chick gut [3]. (g) A simulated soft layer growing in a large simulation domain (77 = 1 and 

9 — Qx — Qy — folds spontaneously into a disordered mixture of Y and T-shaped triple-junction folds and S-shaped zigzag 
folds, reminiscent of the folded cerebral cortex (h) |4T]. 


comparatively ordered near-threshold patterns similar to 
those in Fig. Ef, but the far-from-threshold patterns ap¬ 
pear as disordered mixtures of zigzag-folds and triple¬ 
junctions, reminiscent of the folded cerebral cortex when 
?7 « 1 (Figs. s and|^). The disordered patterns include 
some T-shaped triple-junction folds that lack the rota¬ 
tional symmetry of the Y-shaped ones. T-shaped folds 
do not, however, form a simple periodic pattern and they 
may thus be considered as variants of the ideal Y-shaped 
folds. We note that the disordered patterns resulting 
from unaided folding have slightly higher energies than 
the corresponding periodic zigzag or triple-junction pat¬ 


terns. 


III. SUMMARY 

Our study demonstrates that complex fold patterns 
emerge from even the simplest systems incorporating 
elasticity and differential growth, providing a guide for 
more specific applications incorporating, for example, 
curved geometries, anisotropy, prestressed substrates, 
compressibility and stress relaxation. We show a tran¬ 
sition as a function of modulus ratio from supercritical 
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smooth finite-wavelength buckling followed by sulcifica- 
tion of valleys to direct strongly subcritical sulcification 
at 0.86 < rjs < 1 . 14 , depending on growth anistropy. 
This is accompanied by a transition from contraction to 
elongation of past threshold wavelengths, driving a con¬ 
vergence to spacings within three to four layer thicknesses 
over a range of modulus ratios in both two and three di¬ 
mensions. In equibiaxial growth there is also a beyond 
threshold transition from zigzags to triple-junction folds 
within the typical range of modulus ratios in soft layered 
tissues, underpinning the morphogenetic diversity of me¬ 
chanical folding. 


Appendix A: Details of numerical simulations 

The simulations are based on the explicit finite element 
method with discretizations to constant strain triangle 
( 2 D) or tetrahedron ( 3 D) elements. The elements have 
a quasi-incompressible nodal pressure formulation with a 
bulk modulus K = 10 ^// ( 2 D) or K = 10 ^/i ( 3 D). The 
simulated substrate has a thickness ten times the layer 
thickness. Whether the substrate has a clamped or free 
base has a negligible effect on the results, confirming that 
it provides a good approximation of an infinite substrate. 

In 2 D plane-strain simulations the domain includes a 
half fold with symmetric boundary conditions at sides. 
The domain is discretized to a rectangular mesh of tri¬ 
angles such that the growing layer includes 200 elements 
through its thickness. A fold is initiated by applying 
a small normal force to the surface at the edge of the 
domain. The surface is prevented to cross the vertical 
line that defines the boundary of the domain to prevent 
self-intersection of the fold. The energetically optimal 
wavelength is searched by adjusting the relative width of 
the domain. 

In 3 D simulations the domain has periodic bound¬ 
ary conditions at the lateral boundaries. The domain 
is discretized to tetrahedrons such that the nodes form a 
hexagonal prism mesh. The growing layer includes 10 el¬ 
ements through its thickness. Self-contacts of the surface 
are accounted for by preventing nodes penetrating faces 
at the surface. In the simulations of periodic patterns 
with a unit cell-wide domains point forces are applied 
in selected locations of the domain to initiate formation 
of the patterns. When the dimensions of the domain 
are close to the optimal wavelength of the desired pat¬ 
tern, a stable periodic pattern forms easily. Dimensions 
of the domain are fine-tuned to find the energy minimiz¬ 
ing wavelength. In the simulations of large rectangular 
domains (Fig. |^) the domain has a width 30 times the 
layer thickness and the layer is allowed to fold and re¬ 
lax freely. The surface of the layer in this case has small 
random perturbations that allow folding to initiate. 

The simulation widths for the hysteresis loops in Fig. 
[^were 2 . 50 h for r] = 0 . 5 , 4 . 7 Ih for = I (correspond¬ 
ing to the energy minimizing wavelength ai g = Qs) and 
7 . 54 h for 77 = 4 (corresponding to the wavelength of lin¬ 


ear instability). The hysteresis loops were calculated in 
unloading. The loading curves are deduced from the fact 
we know the system simply becomes unstable at the lower 
of the linear stability threshold or gg. 


Appendix B: Details of linear stability analysis 

The total elastic energy of the layer and substrate is 

/ oo poo pa 

/ / i/uTr (F ■ G-^ ■ G-'^ ■ F^) AzdxAy, 

-00 J —00 J —00 

(Bl) 

with 

I,, r\ - i iv,(^^^&{ 9 x, 9 va/{ 9 xgy))) 0 <z<a 

^<0 

(B 2 ) 

If we restrict attention to x — z plane-strain deformations 
then the above integrand can be written out in terms of 
components of F as 

fciFlv [ghl {FlFFl)Fg-^{FlFFl)) 0 < z < a 
\c2+ I {F^z + Fzz + F^^ -h F^^) z <0 

(B 3 ) 

where the constants ci and C2 do not depend on F, and 
hence do not enter into energy minimization. 

We enforce incompressibility, encoded by Det(F) = I, 
by introducing a Lagrange multiplier pressure field P, 
leading to the effective energy density 

L= |/iTr(P-G“^-G“^-P^)-/iP(Det(P)-I). (B 4 ) 

Minimizing this energy with respect to variations in u 
and P gives the Euler-Lagrange equations 

V • cr = 0 (B 5 ) 

Det(P) = I, (B6) 

where a = dL/d\/u is the stress-tensor, given by 

cr = /i (P . G"^ • G"^ - PP-^) . (B 7 ) 

Energy minimization also requires the free surface to 
be stress free, and the layer and substrate stresses to 
match at their interface: 

o- ■ zU=a = 0 , a- z|^=o+ = o- ■ z|^=o- • (B8) 

We restrict our attention to growth in the x direction, 
9 x = 9 ^ Qy = which we expect buckling in the x — z 

plane so we consider the linear stability of the flat state to 
infinitesimal sinusoidal perturbations with wavenumber 
k, writing the fields as 

P = Po + esm{kx)Pi{z) (B 9 ) 

u = e{cos{kx)d{z)Si + sm{kx) f {z)z). 
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Working in an a; — ^ basis, to first order in e we thus have 


F = I + e 


F-^ = I-e 


—kd{z) sm{kx) cos{kx)d'{z) 
kcos{kx)f{z) sm{kx)f'{z) 


, (BIO) 


—kd{z) sm{kx) kcos{kx)f{z) 
cos{kx)d' (z) sm{kx)f'{z) 


At zeroth order in e we have F = F~^ = I and hence 
piecewise constant stress a = ji {G~^ • G~^ — Pq!)- The 
bulk equations are thus already satisfied, and the bound¬ 
ary conditions require cr • z = 0 throughout, solved by 


Po = • G 


z = 




0 < 2 < a 
z <0. 


(Bll) 


Expanding eqn. |B6| to first order in e reduces it to V • 
u = 0, requiring d{z) = f{z)/k. In the layer, eqn. B5 
expanded to first order requires 


e {fiz) + g^P,{z)) - g^riz) = 0 (B12) 

-g^nz) + efiz) + g^Piiz) = 0. (B13) 


The first of these can be solved algebraically for Pi, giv¬ 
ing 


Pii^) 


k‘^ 



(B14) 


and substituting this into the second gives the governing 
equation for f{z): 


g^fi^) (^) _ (y + 1) ^ ^ o_ 


The same results hold in the substrate, but with g = 1. 
The solution for / in the layer and the substrate is thus 


f( \ _ f< 24 e^^-ha 3 e -haie 0<^<a 

\e^^( 6 i -h zkb2) 2 ; < 0 , 

(B16) 

where the and hi are constants of integration, and 
we have disgarded the unbounded solutions in the sub¬ 
strate. Continuity of z displacement at the interface re¬ 
quires continuity of /2 at 2 : = 0, giving 


hi — CLi F (I 2 F CL 3 F Oja (B17) 

while continuity of x displacement requires continuity of 
f^&t z = Q, giving 

62 = (^2 — — (203 + (J2 + tti)- (B 18 ) 

The first order correction to the free boundary conditions 
at 2 ; = a requires 

«! (5^ + 1) + 02 (y + 1) eW + 2a:ig'^e‘FF’^ 

+ 2 aig^e"FF^ = Q (B 19 ) 

- 201 ^ + 202 ^ 6 ^ - as (y + 1 ) 

+ ai{g^ + l)e"FF'' =Q, (B 20 ) 


and the first order correction to the interfacial boundary 
conditions at 2 ; = 0 require 

Oi (r?y + r? + 2) + 02 {F + ?? - 2) + 203 [F + F 

+ 2aF {g'^V - 1) = 0 (B21) 

2aF{ri + 1 ) - 2a2g‘^{g - 1 ) + 03 [F + F + v) 

+ 04 [F - (/ + l)v)= 0. (B22) 

These four linear equations can be re-expressed as single 
matrix equation of the form M- (ai, a 2 , as, a 4 ) = 0, which 
will only admit non-trivial solutions if Det(M) = 0, the 
condition displayed in eqn. We find the first point of 
instability by finding the wave-number that satisfies this 
condition with the minimum g. Solutions to this eqn are 
shown for a variety of different moduli ratio in fig. la in 
the main text. Each line shows the strain required for in¬ 
stability at each wavelength for a given slab. The actual 
point and wavelength of instability are given by the first 
unstable mode, that is, by the minimum of the curve. 
Inspecting fig. la we see that all layers are unstable to 
a zero wavelength instability ai g = gij = 1.839..., which 
is the Biot strain. We can extract this from the deter¬ 
minant condition by noting that it contains two terms 
which diverge at small wavelength: 

(/ - sy -g^-l) A{g, r,)e“''+7^ + B{g, = 0, 

(B23) 

where 

A={g‘^-F (B24) 

X [g^g^ - y»7(3?7 + 2) - y ( 7 ? + 2f - + 2)) 

B = -(y + i)q/+3y-y + i) (B 25 ) 

X + g^vi^v -‘P)- g‘^(ji- F + {g- F) ■ 

Eor asymptotically small wavelength the left expo¬ 
nential dominates and, since A 7 ^ 0 the con¬ 
dition is only satisfied at the Biot strain ^5 = 

I (yi9 + 3v^+ yi9-3v^+l) = 1.839... where 

g^ — 3^^ — ^^ — 1 = 0. Eurther inspecting fig. la we 
see that, for sufficiently soft layers {g < 0.525...) this 
zero wavelength Biot solution is the first unstable mode, 
whereas for stiffer layers there is a finite wavelength buck¬ 
ling at a lower growth, giving rise to a more conventional 
buckling. We characterize the instability by minimizing 
the curves in fig.la from the main text over wavelength 
to find the first unstable mode for each modulus ratio, 
which gives the wavelength of instability and the growth 
required. These are plotted against modulus ratio in figs. 
Ib-c in the main paper, which again clearly show the 
transition from Biot zero-wavelength behavior to finite 
wavelength buckling at 77 > 0.525.... 

Eor layers marginally stiffer than 775 the threshold for 
instability is just below g^ so we can analyze the insta¬ 
bility by setting g = gb F e and considering the lead¬ 
ing behavior in e. The Biot polynomial then becomes 







— 3^^ — g‘^ — 1 = {pgl — 12gl — 2^^) e so the condition 
for instability becomes 


(656 - 1256 - 256 ) eA{gb,r])e 


cik-\- 


ak-^ 


'>+B{gb,r])e = 0 , 

(B26) 


which is solved by 


_ 2ak 

B{gb,v)e _ 

(656 - 1256 - ^9b)A{gb,g)' 


(B27) 


2)^ + {Vb — 2)7?6 = 0, satisfied if 

_ 56 - 25 ^ + ^/5g^ + 8gf + 2g^ + l + 1 _ „ 

* “ !>S + - !>? + 1 ■ - 

(B28) 

Further analysis indicates that the transition is contin¬ 
uous but logarithmic, in the sense that the unstable 
wavelength vanishes as A ^ 1 / log( 7 ^ — 7 ^ 5 ). 


This small correction to the Biot-threshold for short 
but non-zero wavelengths is positive, indicating the Biot 
point is a minima in fig. la, if 5(^ 5 , 77 ) > 0 and negative, 
indicating the Biot point is a maxima, if 5 (^ 5 , 77 ) < 0. 
The crossover occurs at 5 (^ 5 , 775 ) = 0 which, inspecting 
the form of B, requires g^{gb)^ + dtVbi^Vb - 2 ) - g^gb - 
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